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ABSTRACT 

The algorithm of the ensemble pulsar time based on the optimal Wiener filtration 
method has been constructed. This algorithm allows the separation of the contribu- 
tions to the post-fit pulsar timing residuals of the atomic clock and pulsar itself. Filters 
were designed with the use of the cross- and autocovariance functions of the timing 
residuals. The method has been applied to the timing data of millisecond pulsars PSR 
B1855-I-09 and PSR B1937-I-21 and allowed the filtering out of the atomic scale compo- 
nent from the pulsar data. Direct comparison of the terrestrial time TT(BIPM06) and 
the ensemble pulsar time PTons revealed that fractional instability of TT(BIPM06)- 
PTens is equal to (J:, = (0.8 ± 1.9) • IQ-^^. Based on the a^ statistics of TT(BIPM06)- 
PTens a new limit of the energy density of the gravitational wave background was 
calculated to be equal to ^gh? ^ 3 • 10"^. 

Key words: time - pulsars: individual: PSR B1855-I-09, PSR B1937+21 - methods: 
data analysis 

1 INTRODUCTION 

The discovery of pulsars in 1967 (IHewish et al. 19^8j) showed clearly that their rotational stability 



k> \ allowed them to be used as astronomical clocks. This became even more obvious after discovery 
c^ ' of the millisecond pulsar PSR B1937-I-21 in 1982 (IBacker et al. 19821) . Now a typical accuracy of 



measuring the time of arrivals (TOA) of millisecond pulsar pulses comprises a few microseconds 
and even hundreds of nanoseconds for some pulsars. For the observation interval in the order 
10* seconds this accuracy produces a fractional instability of 10~^^, which is comparable to the 
fractional instability of atomic clocks. Such a high stability cannot but used for time metrology and 
time keeping. 

There are several papers considering applicability of stability of pulsar rotation to time scales. 
The paper (IGuinot 1988P presents principles of the establishment of TT (Terrestrial Time) with 
the main conclusion that one cannot rely on the single atomic standard before authorised confir- 
mation and, for pulsar timing, one should use the most accurate realisations of terrestrial time 
TT(BIPMXX) (Bureau International des Poids et Mesures). The paper of Ilyasov et. al. (I1989P 
describes the principles of a pulsar time scale, definition of "pulsar second" is presented. Guinot & 
Petit (1991) show that, because of the unknown pulsar period and period derivative, rotation of 
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millisecond pulsars is only useful for investigations of time scale stability a posteriori and with long 
data spans. The papers piyasov et al. 1996]), ([Kopeikin 1997[), ([Rodin, Kopeikin &: Ilyasov 1997^, 



( Ilyasov, Kopeikin fc Rodin 1998 ) suggest a binary pulsar time (BPT) scale based on the motion 
of a pulsar in a binary system with theoretical expressions for variations in rotational and binary 
parameters depending on the observational interval. The main conclusion is that BPT at short 
intervals is less stable than the conventional pulsar time scale (PT), but at a longer period of 
observation (10^ 4- 10'^ years) the fractional instability of BPT may be as accurate as 10~^^. The 
paper of Petit & Tavella (119960 presents an algorithm of a group pulsar time scale and some ideas 
about the stability of BPT. The paper of Foster & Backer (ll99Up presents a polynomial approach 
for describing clock & ephemeris errors and the influence of gravitational waves passing through 
the Solar system. 

In this paper the author presents a method of obtaining corrections of the atomic time scale 
relative to PT from pulsar timing observations. The basic idea of the method was published earlier 
in the paper (IRodin 2006^ . 

In Sect. [2], the main formulae of pulsar timing are described. Sect. [3] contains theoretical al- 
gorithm of Wiener filtering. Sect. [5] presents the results of numerical simulation, i.e. recovery 
of harmonic signal from noisy data by Wiener optimal filter and weighted average algorithm. 
The latter one is used similarly to the paper of Petit & Tavella (119961) . Sect. [6] describes an 
application of the algorithm to real timing data of pulsars PSR B1855+09 and PSR B1937+21 



dKaspi, Taylor fc Ryba 19941). 



2 PULSAR TIMING 

The algorithm of the pulsar timing is widely described in the literature ([Backer fc Hellings, 1986[) 



([Doroshenko fc Kopeikin 1990P, ([Edwards, Hobbs fc Manchester 2d06|). Two basic equations are 



presented below. Expression for the pulsar rotational phase (/)(t) can be written in the following 
form: 

0(t) = 0o + z/t + -z>t2 + £(t), (1) 

where t is the barycentric time, 0o is the initial phase at epoch t = 0, z/ and z> are the pulsar spin 
frequency and its derivative respectively at epoch t = 0, and e{t) is the phase variations (timing 
noise). Based on the eq. ([T]) pulsar rotational parameters v and z> can be determined. 

The relationship between the time of arrival of the same pulse to the Solar system barycentre 
t and to observer site t can be described by the following equation (Murray, 1986) 

c(t - t) = -(k ■ b) + — [k X hf + Atrci + At^M, (2) 
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where k is the barycentric unit vector directed to the pulsar, b is radius- vector of the radio telescope, 
i? is a distance to the pulsar, Atrei is the gravitational delay caused by the space-time curvature, 
At^M is the plasma delay. The pulsar coordinates, proper motion and distance are obtained from the 
eq. ([2]) by fitting procedure that includes adjustment of above-mentioned parameters for minimising 
the weighted sum of squared differences between 0(t) and the nearest integer. 

3 FILTERING TECHNIQUE 

Let us consider n uniform measurements of a random value (post-fit timing residuals) r = 
(ri,r2, . . . ,rn). r is a sum of two uncorrelated values r = s + e, where s is a random signal to 
be estimated and associated with clock contribution, e is random error associated with fluctuations 
of pulsar rotation. Both values s and e should be related to the ideal time scale since pulsars in the 
sky "do not know" about time scales used for their timing. The problem of Wiener filtration is con- 
cluded in estimation of the signal s if measurements r and covariances ([3]) are given (1 Wiener 1949t 
IGubanov 1997^ . For r, s and e their covariance functions can be written as follows 



{r,r) = {ri,rj) = {si,Sj) + {ei,ej) 



cov 

cov(s,s) = (si,s,), fi 7 = 1 2 n) (3) 

coy{s,r) = {si,rj) = {ri,Sj) = {si,Sj) , ^ '-^ ''■■■'' 
cov(e,£:) = {suSj) . 

denominates the ensemble average. 

The optimal Wiener estimation of the signal s and a posteriori estimation of its covariance 

function D^s are expressed by formulae (IWiener 1949[ IGubanov 19970 

S = QsrQ,T> = Q,,Q;> = QssiQss + Qee)"^r (4) 

and 

where covariance matrices Q^r, Qsr, Qrs, Qss are constructed as Toeplitz matrices from the corre- 
sponding covariances. In this paper we assume that processes s and e are stationary in a weak sense 
(stationary are the first and second moments). Since quadratic fit eliminates the non-stationary part 



of a random process (Kopeikin 1999), their covariance functions depend on the difference of the 



time moments tj — tj. 

Since it is impossible to perform pulsar timing observations without a reference clock, for sepa- 
ration of the covariances, (sj, Sj) and {ei,ej), it is necessary to observe at least two pulsars relative 
to the same time scale. In this case, combining the pulsar TOA residuals and accepting that cross- 
covariances (^Ei, ^Ej) = {^Ei,^Ej) = produces (hereafter upper-left indices run over pulsars under 
consideration) 
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{si, Sj) = ((|Vi + Vi, Vj + Vj^ - (|Vi - Vj, Vj - Vj^) /4, or (si, Sj) = (|Vj, Vj^ . (6) 

If M pulsars are used for construction of the pulsar time scale then one has M{M — l)/2 signal 
covariance estimates (s,, Sj) = l^Ti, Wj\ , (A;, / = 1, 2, . . . , M). 

In formula (jl]) the matrix Q~^^ serves as the whitening filter. The matrix Q^s forms the signal 
from the whitened data. 

The ensemble signal (pulsar time scale) is expressed as follows 

J/(M-1) 

2 2 M 

Sens = TTTTF^TTT ^ m^ss ' 2_^ W (4^^ ■ T, (7j 

where % is the relative weight of the ith pulsar, *w = n/crf., Uj is the root-mean-square of whitened 
data *Q,~^ ■ *r, and k is the constant serving to satisfy Y^i^w = 1. The first multiplier in eq.© is 
the average cross-covariance function, the second multiplier is the weighted sum of the whitened 
data. 

For calculation of the auto- and cross-covariances, the following algorithm was used: the initial 
time series ^rt were Fourier-transformed 

1 " 
'x{u) = — E 'n/ite^^*, (A; = 1, 2, . . . , M), (8) 



weights ht are the 0th order discrete prolate spheroidal sequences (IPercival 199ip which are used 



for optimisation of broad-band bias control. These can be calculated to a very good approximation 



using the following formula (IPercival 199ip 



h[7:W{n-l)^ll-{^-iy 

ht = a- 



(9) 



loinWin-l)) 

where Cq is the scaling constant used to force the convention J2hf = 1, Iq is the modified Bessel 
function of the 1st kind and 0th order, and the parameter W affects the magnitude of the side-lobes 
in the spectral estimates (usually W = 1 -r 4). In this paper W = 1 was used. 

Power spectrum {k = I) and cross-spectrum {k ^ I) were calculated by the formula 

1 



'''X(uj) - — 
^ ' 27r 



x{LUyx*{uj) , (10) 



where (■)* denominates complex conjugation. 

Lastly, auto- {k = I) and cross-covariance {k ^ I) were calculated using the following formula 



covf'^r/r) 



Y^''X{u;)e-'^\ {k,l = 1,2,..., M). (11) 



UJ = 1 
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Figure 1. The accuracy of signal estimation based on the methods of weighted average (dashed line) and Wiener filter (solid line) as 
dependent on the number of pulsars (left panels) and length of the data (right panels). For the calculation shown in the left panels 256 
points of data were taken, for calculations shown in the right panels five pulsars were used. Different types of noise were generated: (a), 
(b) - white phase noise, (c), (d) - white noise in frequency, (e), (f) - random walk noise in frequency. Data in the panels (d) and (f) were 
scaled accordingly for fitting within in all the panels. 



4 COMPUTER SIMULATION 

To evaluate performance of the Wiener filtering method as compared to the weighted average 
method, we have applied it to simulated time sequences corresponding to harmonic signal with 
additive white and red (correlated) noise. Simulated time series were generated with the help of 
random generator built in the Mathematica software which had the preset (normal) distribution for 
different numbers of pulsars. A maximum of 50 pulsars were used (limited by acceptable computing 
time). The harmonic signal to be estimated was as follows: Asin{t/P), A = 1, P = 10, t = 
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1,2,..., 256. The additive Gaussian white noise had zero mean and different variance. For example, 
in simulation for 50 pulsars, the variance was in the range of a^ = 1, 2, . . . , 50. The correlated noise 
^2, n^ with the power spectra 1/ P and 1//^ was generated as a single or twice repeated cumulative 
sum of the white noise hq: 

j j 

'n'2j = Y.^oi^ '^'ij = J2'^2i, (j = l,2,...,n). (12) 

The second order polynomial trend was subtracted from the generated time series n2 and n^. 
The weights of the individual time series were taken inversely proportional to cr^, where o"^ is the 



fractional instability (Taylor 1991). Quality of the two methods was compared by calculating the 
root mean square of the difference between original and recovered signals. 

Fig. [T] shows the results of computer simulation described above. Quality of these two methods 
of signal estimation is clearly visible. It is important to note that signal estimation accuracy in the 
case of Wiener filter (solid line) is better in all cases. The most significant advantage of the Wiener 
filter over the weighted average method is seen in the case of the white noise (fig. 1(a), 1(b)). For 
the correlated noise with the power spectrum 1//^ (fig. 1(c), 1(d)) the advantage is still clear. For 
the red noise with the power spectrum 1//^ both methods show similar results (fig. 1(e), 1(f)). 
Noteworthy is dependence of estimation quality on the observation interval for the correlated noise 
(fig. 1(d), 1(f)): as the observation interval increases, the estimation accuracy grows. Physically such 
a behaviour corresponds to appearance of more and more strong variations of the correlated noise 
with time, which deteriorate the signal estimation quality. Influence of the form of the correlated 
noise and length of the observation interval on the variances of the pulsar timing parameters are 



described in detail in (Kopeikin 1997) 



5 RESULTS 



To evaluate the performance of the proposed Wiener filter method, timing data of pulsars PSR 
B1855+09 and B1937+21 ([Kaspi, Taylor &: Ryba 1994^ were used. For the sake of simplicity of 



the subsequent matrix computations, unevenly spaced data were transformed into uniform ones 
with a sampling interval of 10 days by means of linear interpolation. Such a procedure perturbs a 
high-frequency component of the data while leaving a low- frequency component, of most interest 
to us, unchanged. The sampling interval of 10 days was chosen to preserve approximately the same 
volume of data. 

The common part of the residuals for both pulsars (251 TO As) has been taken within the 
interval MJD = 46450 -^ 48950. Since choosing the common part of the time series changes the 
mean and slope, the residuals have been quadratically refitted for consistency with the classical 
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Figure 2. The barycentric timing residuals of pulsars PSR B1855+09 (thin line) and PSR B1937+21 (solid line) before (a) and after 
(b) uniform sampling. 



timing fit. The pulsar post-fit timing residuals before and after processing described above are 
shown in Fig. [2l 

According to ( |Kaspi, Taylor fc Ryba 19941 ), the timing data of PSRs B1855+09 and B1937+21 
are in UTC time scale. UTC (Universal Coordinated Time) is the international atomic time scale 
that serves as the basis for timekeeping for most of the world. UTC runs at the same frequency 
as TAl (International Atomic Time). It differs from TAI by an integral number of seconds. This 
difference increases when so-called leap seconds occur. The purpose of adding leap seconds is to keep 
atomic time (UTC) within ±0.9 s of an time scale called UTl, which is based on the rotational rate 
of the Earth. Local realisations of UTC exist at the national time laboratories. These laboratories 
participate in the calculation of the international time scales by sending their clock data to the 
BIPM. The difference between UTC (computed by BIPM) and any other timing centre's UTC 
only becomes known after computation and dissemination of UTC, which occurs about two weeks 
later. This difference is presently limited mainly by the long-term frequency instability of UTC 
flAudoin fc Guinot 20011 ). 

The signal we need to estimate is the difference UTC - PT. Fig. [3] shows the signal estimates 
(thin line) based on the timing residuals of pulsars PSR B1855+09 and PSR B1937+21 and cal- 
culated with use of eq. (jl]). The combined signal (ensemble pulsar time scale, eq. ([7])) is shown in 
fig. m and displays behaviour similar to the difference UTC - TT(BIPM06) (correlation p = 0.75). 

All three signals UTC - PT1855, UTC - PT1937 and UTC - PTens were smoothed by use 
of the following method: to decrease the Gibbs phenomenon (signal oscillations) near the ends 
of smoothing interval, the series under consideration were backward and forward forecasted by 
p = Integer Part[n/2] lags {n = 251 is length of time series) with use of the Burg's (also re- 



ferred to as the maximum entropy method) autoregression algorithm of order p ( Burg, 1975 
ITerebizh 19920 . A new time series of double length were smoothed by use of the low-pass Kaiser fil- 
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Figure 3. Differences in UTC - PT1855 (a) and UTC - PT1937 (b) (thin line) for the interval MJD = 46450 -^ 48950 estimated using 
the optimal filtering method (eq. (4)). The smoothing with Kaiser filter is shown by solid line. The dashed line displays the difference 
UTC - TT{BIPM06). 
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Figure 4. Combined clock variations of UTC — PTcns for the interval MJD = 46450 -\- 48950 estimated using the optimal filtering 
method from the timing residuals of pulsars PSR B1855+09 and PSR B1937+21 (solid line) and difference UTC - TT(BIPM06) (dashed 
line). The thin line indicates the difference TT - PTons 



ter fIGold fc Rader 1973HKaiser 19741) with the bandwidth of /max/32, where /^ax = 2/ At, At = 10 
days is the samphng interval. The choice of the bandwidth was defined by visual comparison with 
the UTC - TT(BIPM06) line. The final time series were obtained by dropping the forward and 
backward forecasted sections of the smoothed series. 

The accuracy of the obtained reahsations of the pulsar time UTC - PT1855 and UTC - PT1937 
was derived from the diagonal elements of the covariance matrix defined by eq. ([5]). The root mean 
square value of UTC - PT1855 and UTC - PT1937 is equal to 0.44 yus and 0.67 /is respectively. The 
accuracy of the smoothed signals was estimated as 0.44/vT6 = 0.11 yus and 0.67/vT6 = 0.17 /is. 
Finally, for overall accuracy a conservative estimate 0.17 /xs was accepted. 
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6 DISCUSSION 

The stability of a time scale is characterised by so-called Allan variance numerically expressed as a 
second-order difference of the clock phase variations. Since timing analysis usually includes deter- 
mination of the pulsar spin parameters up to at least the first derivative of the rotational frequency, 
it is equivalent to excluding the second order derivative from pulsar TOA residuals and therefore 
there is no sense in the Allan variance. For this reason, for calculation of the fractional instability 



of a pulsar as a clock, another statistic cr^ has been proposed (Taylor 1991). A detailed numerical 



algorithm for calculation of o"^ has been described in the paper ([Matsakis, Taylor &: Eubanks 1997|) 



In this work, an idea has been proven that different realisations of pulsar time scales must 



have comparable stability between each other (Lyne fc Graham-Smith, 1998) and should be not 



worse than available terrestrial time scale at the same interval. For this purpose statistic a^ of two 
realisations of PT, UTC - PT1855 and UTC - PT1937, were compared. 

Fig. E] presents the fractional instability of the differences PT1937 - PT1855 (dashed line) and TT 
- PTens (sohd line). At 7 years time interval a, = (0.8 ± 1.9) ■ 10"^^ and a, = (1.6 ± 2.9) • 10"^^ for 
TT - PTens and PT1937 - PT1855 respectively. One can see that instability of the two differences lays 
within error bar intervals. The fractional instability of TT relative to PTens and PT1937 relative to 
PTi855 is almost one order of magnitude better than individual fractional instability of the pulsars 
PSR B1855+09 and PSR B1937+21. 

As an example of astrophysical application of the fractional instability values obtained in this 
work, one could consider estimation of the upper limit of the energy density of the stochastic 
gravitational wave background (Kaspi, Taylor & Ryba 1994). For this purpose theoretical lines of a^ 



in the case when the gravitational wave background with the fractional energy density Qgh"^ = 10~^ 
and 10"^'' begins to dominate are plotted in the lower right hand side corner of fig. O One can see 
that a^ of TT - PTens crosses the line Qgh'^ = 10~^ and approaches Qgh'^ = 10^^°. The upper limit 
of Qgh"^ based on the two sigma uncertainty (95% confidence level) of the ensemble az is equal to 
- 3 ■ 10-9. 

Noteworthy, that since PSRs 1855+09 and 1937+21 are relatively close to each other in the sky 
(angular separation is 15.5°) and hence their variations of the rotational phase contain a correlated 
contribution caused by the stochastic gravitational wave background (Hellings & Downs 1983), they 
form a good pair for estimation of the upper limit of Qgh"^. 

Currently, the accuracy of the filtering method without contribution of the uncertainty of the 
TT algorithm is estimated at a level of 0.17 fis. So, the uncertainty in PTens may, in principle, 
reach the level of a few ten nanoseconds if it were to be used for all high-stable millisecond pulsars. 
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Figure 5. The fractional instability (Jz based on the difference PT1937-PT1855 (dashed line) and cfz of the difference TT — PTens (solid 
line). Theoretical cfz in the case Qgh^ = lO^'-* and 10^^" are shown in the lower right-hand corner of the plot. 

As computer simulations show, for the highest advantage while applying the Wiener optimal filters 
one should use the pulsars that show no correlated noise in their post-fit timing residuals. 



7 CONCLUSIONS 

An algorithm of forming of ensemble pulsar time scale based on the method of the optimal Wiener 
filtering is presented. The basic idea of the algorithm consists in the use of optimal filter for removal 
of additive noise from the timing data before construction of the weighted average ensemble time 
scale. 

Such a filtering approach offers an advantage over the weighted average algorithm since it 
utilises additional statistical information about common signal in the form of its covariance function 
or power spectrum. Since timing data are always available relative to a definite time scale, for 
separation of the pulsar and clock contributions one need to use observations from a few pulsars 
(minimum two) relative to the same time scale. Such approach allows estimation of the signal 
covariance function (power spectrum) by averaging all cross-covariance functions or power cross- 
spectra of the original data. 

The availability of two scale differences UTC - TT and UTC - PT has resulted in the long 
awaited possibility of comparing the ultimate terrestrial time scale TT and extraterrestrial ensemble 
pulsar time scale PT of comparable accuracy. The fractional instability of the terrestrial time scale 
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TT relative to PT and their high correlation have demonstrated that PT scale can be successfully 
used for monitoring the long-term variations of TT. 
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